{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discussion 3: Probabilistic View of Linear Regression\n",
    "\n",
    "[CS189/289A - Fall 2020]\n",
    "\n",
    "## Introduction\n",
    "\n",
    "Welcome! This discussion material will cover parts of Problem 2 and 4 in HW3. It aims to give you a better understanding of an interpretation of regression as a probabilistic model and should get you started on the homework problems. We will first start with this Jupyter notebook and play around with multivariate Gaussian distribution (2D). We will consider a linear regression problem from a probabilistic view and examine the likelihood function and the posterior function.\n",
    "\n",
    "## (a) Isocontour of Gaussian PDF\n",
    "\n",
    "We will warm-up by getting all of you more familiar with multivariate Gaussian distribution by observing the contour of its density in 2D. Run the following code and use adjust the following parameters: mean of x and y coordinates (`mu_x` and `mu_y`), and its covariance matrix, $\\Sigma \\in \\mathbb{R}^{2 \\times 2}$, given by\n",
    "\n",
    "$$\n",
    "\\Sigma = \\begin{bmatrix}\n",
    "\\Sigma_{xx} & \\Sigma_{xy} \\\\\n",
    "\\Sigma_{xy} & \\Sigma_{yy}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "For some values of $\\Sigma$, you may get a matrix that is no longer positive semidefinite. If that happens, move $\\Sigma_{xy}$ closer to zero or just try a different set of values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "from scipy.integrate import simps\n",
    "from ipywidgets import interactive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0e16e8bb25f14d3eb445a7e14ec9554e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, description='mu_x', max=4, min=-4), IntSlider(value=0, description='m…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_contour(mu_x, mu_y, s_xx, s_yy, s_xy):\n",
    "    # Generate grid of points at which to evaluate pdf\n",
    "    x = np.linspace(-5, 5, 500)\n",
    "    y = np.linspace(-5, 5, 500)\n",
    "    X,Y = np.meshgrid(x, y)\n",
    "    pos = np.array([Y, X]).T\n",
    "    S = [[s_xx, s_xy], [s_xy, s_yy]]\n",
    "    rv = scipy.stats.multivariate_normal([mu_x, mu_y], S)\n",
    "    Z = rv.pdf(pos)\n",
    "    \n",
    "    e, v = np.linalg.eig(S)\n",
    "    v = v.T * np.sqrt(e.reshape(2, 1))\n",
    "    plt.figure(figsize=(6, 5))\n",
    "    plt.contourf(X, Y, Z)\n",
    "    plt.arrow(mu_x, mu_y, v[0, 0], v[0, 1], width=0.1, color='r')\n",
    "    plt.arrow(mu_x, mu_y, v[1, 0], v[1, 1], width=0.1, color='r')\n",
    "    plt.colorbar()\n",
    "    plt.xlim(-5, 5)\n",
    "    plt.ylim(-5, 5)\n",
    "    plt.xlabel('x')\n",
    "    plt.ylabel('y')\n",
    "    plt.show()\n",
    "\n",
    "interactive_plot = interactive(plot_contour, mu_x=(-4, 4), mu_y=(-4, 4), \n",
    "                               s_xx=(0.1, 5), s_yy=(0.1, 5), s_xy=(-5, 5))\n",
    "output = interactive_plot.children[-1]\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer the following questions:\n",
    "\n",
    "- __What happens when you change the mean?__: Changing the mean translates the Gaussian PDF.\n",
    "- __What happens when you change the diagonal entries of the covariance matrix? (also notice the color bar on the right)__: Modifying the covariance matrix linearly transforms the PDF. Changing the diagonal entries while keeping the off-diagonal zero scales the PDF along x and y axes.\n",
    "- __What happens when you change the off-diagonal entries?__: Doing so rotates the contour of the PDF. \n",
    "\n",
    "The red arrows show the two orthogonal directions with the largest and the smallest variance. They can be seen as the transformed coordinates! If you have data points drawn from a standard Gaussian distribution (i.e. $\\Sigma = I$), they can be linearly transformed as if they are drawn from any Gaussian distribution with covariance matrix $\\Sigma$. The transformation matrix is given by \n",
    "\n",
    "$$T = VS^{1/2}$$\n",
    "\n",
    "where $V$ is the eigenvectors of $\\Sigma$, and $S$ is the diagonal matrix of eigenvalues of $\\Sigma$. In other words, SVD of the covariance matrix $\\Sigma$ is given by\n",
    "\n",
    "$$\\Sigma = VSV^\\top = TT^\\top$$\n",
    "\n",
    "The red arrows are then generated by transforming the standard unit vectors, $[0, 1], [1, 0]$, by the transformation matrix $T$.\n",
    "\n",
    "__Why does this make sense mathematically?__ Consider PDF of an standard Gaussian: $f(x) = \\frac{1}{\\sqrt{2\\pi}} \\exp\\left\\{-\\frac{1}{2}\\|x\\|_2^2\\right\\}$. We want to find PDF of the transformed distribution $g(z)$ where $z = Tx$. More generally, for any PDF $f(x)$, $g(z)$ is given by\n",
    "\n",
    "$$g(z) = f(x(z)) \\left| \\det\\left( \\frac{dx}{dz} \\right)\\right|$$\n",
    "\n",
    "Now see what $g(z)$ looks like for our Gaussian distribution.\n",
    "\\begin{align}\n",
    "g(z) &= \\frac{1}{\\sqrt{2\\pi}} \\left| \\det{T^{-1}} \\right| \\cdot \\exp\\left\\{-\\frac{1}{2}\\|T^{-1}z\\|_2^2\\right\\} \\\\\n",
    "&= \\frac{1}{\\sqrt{2\\pi}} \\left| \\frac{1}{\\det{T}} \\right| \\cdot \\exp\\left\\{-\\frac{1}{2} z^\\top (T^{-1})^\\top T^{-1} z \\right\\} \\\\\n",
    "&= \\frac{1}{\\sqrt{2\\pi}} \\left| \\frac{1}{\\sqrt{\\det{\\Sigma}}} \\right| \\cdot \\exp\\left\\{-\\frac{1}{2} z^\\top (S^{-1/2}V^\\top)^\\top S^{-1/2}V^\\top z \\right\\} \\\\\n",
    "&= \\frac{1}{\\sqrt{2\\pi(\\det{\\Sigma})}} \\cdot \\exp\\left\\{-\\frac{1}{2} z^\\top V S^{-1}V^\\top z \\right\\} \\\\\n",
    "&= \\frac{1}{\\sqrt{2\\pi(\\det{\\Sigma})}} \\cdot \\exp\\left\\{-\\frac{1}{2} z^\\top \\Sigma^{-1} z \\right\\} \n",
    "\\end{align}\n",
    "\n",
    "This shows that transforming $x$ by $T$ is equivalent to sampling $z = Tx$ from a new Gaussian distribution with covariance matrix $\\Sigma$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (b) Linear Regression (Likelihood)\n",
    "\n",
    "Given the following variables: a data point $x \\in \\mathbb{R}^d$ and its label $y^* \\in \\mathbb{R}$. $x$ and $y^*$ are related by the following linear function.\n",
    "\n",
    "$$y^* = x^\\top w^*$$\n",
    "\n",
    "However, in this case, we can only observe a noisy version of $y^*$:\n",
    "\n",
    "$$y = y^* + \\epsilon = x^\\top w^* + \\epsilon$$\n",
    "\n",
    "where $\\epsilon \\sim \\mathcal{N}(0, \\sigma^2)$. In other words, we are given a training set of $n$ samples: $\\{(x_1,y_1),\\dots,(x_n,y_n)\\}$. This is the same setting as in Problem 2 of HW3 where we will later show the probabilistic interpretation of least square and ridge regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set style for prettier plots\n",
    "plt.style.use('seaborn-darkgrid')\n",
    "\n",
    "def generate_data(n, d, sigma):\n",
    "    np.random.seed(1)\n",
    "    w = 1.\n",
    "    X = np.random.rand(n) * 10 - 5\n",
    "    y = X.dot(w) + np.random.randn(n) * sigma\n",
    "    return y, X, w\n",
    "    \n",
    "num_samples = 200\n",
    "d = 1\n",
    "sigma = 1\n",
    "y, X, w = generate_data(num_samples, d, sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEBCAYAAACdctWRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAApZ0lEQVR4nO3de3BV5fkv8O++JHsDCYk6kED7GxCOQYYO2DKCDg1UbMrFoj0WBVRmGFupYquWFHRAIl5AOJKe6UwPUGx1mBZbhwP8CCOh5dIKpjPoWKGHTgy3gZmaEFPFXCA72TvZ54+44r6877rsvdbaa639/fylyc667JBnvft5n/d5ffF4PA4iIvIkf64vgIiIrMMgT0TkYQzyREQexiBPRORhDPJERB7GIE9E5GHBXF9Aora2zlxfQkaKikLo6urJ9WXYivecH3jP7jBiRLH0exzJmyAYDOT6EmzHe84PvGf3Y5AnIvIwBnkiIg9jkCci8jAGeSIiD2OQJyLyMEeVUBIReUV9Yyu2nriE1s4elBWHsKJyLOZNLLP9OhjkiYhMVt/Yio1/OYdIrB8AcKWzBxv/cg4AbA/0TNcQEZls64lLgwFeEYn1Y+uJS7ZfC0fyREQ66U3BtHaKV8zKvm4ljuSJiHRQUjBXOnsQx1cpmPrG1rTXlhWHhMeQfd1KDPJERDoYScGsqByLcDA5vIaDfqyoHGvhFYrZmq75zW9+g2PHjiEajWLJkiV44IEH7Dw9EZEuorSMkRSMksLJq+qakydP4qOPPsIf//hHdHd344033rDr1EREutWdbhZWxgwPB9EeiaW9XpaCSQ30yojf7kBvW5B/7733UFFRgSeffBJdXV1YvXq1XacmItKt9vBZYVoG8Xjaa9VSME4po7QtyF+9ehXNzc3Yvn07/v3vf+OJJ57AoUOH4PP5Bl9TVBRyZZvPQMCP0tKhub4MW/Ge80M+3nNLe0T49UhfepC//1tfw5I7bxa+fnvDZeHDYnvDZenPWMG2IF9aWopx48ahsLAQ48aNQygUwueff46bbrpp8DVua9SvKC0dii++uJ7ry7AV7zk/5OM9jyoJo1kS6FMda/wUP68UB2zZw6KlPWL6e+qITUOmTp2KEydOIB6Po7W1Fd3d3SgtLbXr9ERESeobW7Fgx0lMqz2OBTtODpZCVldVpFXGyKjVvTuljNK2kfxdd92FDz74AAsXLkQ8HkdNTQ0CAfelZojI/dTy5UvuvBnXrvckVcZ0R/uEk64+38CxRDn2FZVjk84B5KaM0hePC2YTcsSte7zm40da3nN+8Oo9L9hxElcEo/Dy4hBOrL4r7Z5THwqJwkE/1nzvFmGgt6tJmVq6hm0NiCjvGG07oATm9fVN6E8ZFisLoowEbzs7VDLIE1HeKSsOCUfyavnyeRPL8MLBJuH3RA8HWUro9CfteOdfn9pWWsm2BkSUdzJtO6BnMlWZ0K052CQsodz3zyu2dqhkkCeivDNvYhnWfO8WlBeH4MNALl6WV0+k9XBIbGImk5ruUVjVoZLpGiLKS/MmlhlOj2j1pBE1MUvl94kDvVWllQzyROQJdk1mqj0ctEbj4aAf90wamZSTV75uVWklgzwRuZ5T+sTIJnSBgZSQ8uCZ8rUSVtcQEeml1uvdziAvWgAlkkmqKFMM8kTketlst5ea5lk1ZwJmjinN6DoSc/OpI/rETxfKaziSJyLSIZO6d0Cc5lm7/wzWVMkrbbRy/8ooXbSqNhLrx5aj59HbF2edPBE5l6y5V65kWvcuTPNE5TXrRvZ5lX2K6OjpY508ETmXkUBn9XUoD5qtJy7hnkkjDde9ywLxlc4e4cPLyD6vRksiWSdPRI7ghElOUZrlnX99qiuwJ1KrhhGlUYzk/mVdKENBv6FtBLPFkTwRGZLNJKdZjIyo1YjSPGrHlAXi4lAgLX0lW1VbPXt8RqmlTHEkT0SGZDrJaSY9Dxo9i6PUqmFExxSNzoM+oDvaj46egdelVtGk0lo1azYGeSIyxAmbYWg9aIwsjlKrhkk8ZuLPam0oEon1o/bYBfTE+qXXYFdqi+kaIjIk0+ZeZpKlWbqjfYMjeKPpHL0VOvMmluHA8ul4v3omDiyfjg5Bfh0A2iMxW6toZDiSJyLD7ByJys4PAFuOnkdHT9/g19sjMdUVp2rzBsoxtzdcRkt7RHcaRW3y1ug1WIFBnogcR28+feuJS0lBHhgYLRvt9Jh4vlElYbw4f4Luh5gsfVUY8KVdm9o1WIVBnsjh7NwqzgmM5NNlo+L++ECg1TNvkHq+5vaIoRWosolUADmfuwByEOQ/++wz3H///XjjjTcwfvx4u09P5CpO6a5op9pjF3TX4ctSJUrHRz0PR1n+fn39wFZ/egO97HW5fkDbGuSj0ShqamoQDoftPC2Razlh4ZGd6htbhQuFAGMLjpRgquc9Uvs0kO0DVc81WP1JzdYgv3nzZixevBg7duyw87REruWEhUd2Uqs8EeWyzag5V5s4TayGUTtHpoHajk9qtgX5vXv34sYbb0RlZaU0yBcVhRAMBuy6JNMEAn6Ulg7N9WXYivdsj1ElYTS3R4Rft+Na7L5ntYfXqjkThNey5M6bseTOmzM+56o5E7B2/xlEouKKnCudPdh4+Nzg95X/HzY0hHunjEbd6WbV76vZ3nBZ+Elte8PlrO4pkS8ej0u2lTXXww8/DJ/PB5/Ph8bGRowdOxbbtm3DiBEjBl/T1tZpx6WYrrR0KL744nquL8NWvGd7pI70gIF0hF116Xbfs2xB0vBQAEd/OsPw8fSOsOsbW7G+vkm6ybZIeXEIB5ZPl16z8n0102qPQ3RKH4D3q2fqvpYRI4ql37NtJL9r167B/166dCnWr1+fFOCJKJ3dS+BzTZZj/8Xd/8PwsYyuegXSq2HUKIE9m5SaHS0iWEJJ5HC5XnhkJzMfampVMy8cbEo79ryJZTj9STv2nL6i+xz1ja1ZBWo7WkTkJMj//ve/z8VpicgFzHqoqVXNAOKRfcPFq4bOsfXEpawCtR2f1DiSJyJP0tNuILUc1WjVUmtnT9aB2upPagzyRGQ6J6zSFY2wRRIDu9E+NEpKxskpNXahJCJTOWV7wNRumX6f+HU+HwavbUXlWIQL0sPi7f813NaNPszEkTwRmcqsVbpmfBpIHGGLylGB9JWtw4aG8Nqfm9LOK7seJ3xqUcMgT0SmMmOVbn1jK16qb0IsYZL0JQO9ZETUqmcSH0L3ThmNmWNKhT+fem439BZiuoaITCUrHTRS+73l6PnBAK+IxQe+nqn6xla8869Ppd/PpFWEWXvNWolBnohMpXeHJTWiPuxqX9dDFJATZbIAyQ29hZiuISJTOXWVrlrgzXQS1QmbmmthkCci02VbUlgSDgpbDpeEMw9ZsoDs9yHjXkBO2NRcC4M8EelmRSWJ6JjVs8fj5UNnEU3oGFbg96F6duYbDckCcjbN3pz6qSWRbV0o9WAXSvfgPbtPJgE68Z71dsQ0ch61YwLi4JnNg0bPz7rx9+yILpRElDtmlPrpqX/fdORsUomi1nnUjnlg+XRXliw6DatriPKAGaV+WpUk9Y2tqjXomRwzVTb34ZSVuHZjkCfKA2aU+mnVv6sFWtl5jNbUZ3MfbqhptwKDPFEeMGOBklb9u1qglZ3HaE19Nvfhhpp2KzDIE+UBrWBa39iKBTtOYlrtcSzYcVKYwkht+FVeHEqadFULtLKgrXVMo/ehxowHnRtx4pUoD6iV+okmM2sONuH0J+3YtPC2tOOoBWBRA7Ahgq6Oeo9p5D60uKGm3QosoTSBG0uussV79gatDaxrF07GzDGlhjbE3nL0fFr7ASP16FZ2dczHEkoGeRO48R9FtnjP7idrvZtodEkYj88YY2gR0YIdJ4UrS8uLQziwfLrha8p2wZJRbvw9qwV55uSJ8pRWwy4AaGmPGK5KYQWMs9iWk49Go1izZg0++eQT9Pb24oknnsDdd99t1+mJ6EtKykLPNnejSsJoaY8Iv6dWFplp0658rYCxkm0j+bq6OpSWluKtt97Cb3/7W7z88st2nZrIVnoqVXIlcUGQlnDQj+qqCsNVKayAcRbbRvJz587FnDlzAADxeByBQMCuUxPZxo5l99lMTOpJ0QAD+fMVlWNx75TRuHa9x3BVSmHAB6WJZEk4iOrZ41kBkyO2Bflhw4YBALq6uvDUU0/hmWeeSXtNUVEIwaD7gn8g4Edp6dBcX4ateM9i2xsuC3PK2xsuY8mdN2d9DXWnm7Hx8DlEosnljh//5zpeXDBJ8+fV0h6jS8KorqrAvVNGD34tEPBjyZ03Y9jQEGoPn0VLewSjBK+TXR8AtEdi+Pg/17FEx78XI+eyitf+bdtaXdPS0oInn3wSDz30EBYuXJj2fVbXuAfv+SuJI2vZH5MPwPvVM1WPr2eELqtcAYCX5k/QHC3Lfr4kHMSQgkDauY3+nrO9Pidw479tR1TX/Oc//8Gjjz6KVatWCQM8kRulNr2S0cop622epTYS11OBIsqXF/h9uNYTM6Vxl1qunxUyuWFbkN++fTs6OjqwdetWLF26FEuXLkUkIp61J3ILPTlurZyysiBJT+mg2sNCTwWKqI3AkAJ/2qbZmZQtaj0UWCGTG7bl5J9//nk8//zzdp2OyBZqgcsHSNMuessYU4+/onIsag42CV+r59NCYjroxS/TJ9Nqj+s6txathwIrZHKDvWuIsiCrCVdb3alnpWni8RPNm1iG05+0p/Vt1/NpQVb1o6euXc98gRUbZVP2uOKVKAuZ1ITrLWOUHee571bgpfkTdHdulJ1TScmI7gEYeBDM2vI3bDpyVtd8gWykns1G2ZQ9juSJsqCnK2LqKFjPQiStwCjr3CgbcautJE28h9Rra26PqO72lHgNsi6URYXuK4v2EgZ5oiyptcoVpUm0ZNqQK5uUjHIPaiWQqVIfHMr11h67gHZlJRSAjp4+7sOaQ0zXEFlIb2pGoSf1YuRcaikZUToo2+0A500sw5CC9JE7m4zlDkfy5FpW9h03i1rQLC8OZXTt2aZk1M6pN52kNu/AJmPOwiBPrmRHjxgzZFJ9o0a2i1PNwSYMKfCjO5r+qSE1JaNG1jvmnkkj0XDxqq6Hkuyei0PMzecCgzy5klpqwklB3uyGW2rpn+5oPwI+oC9hYZPRc6WO+Ed9uWmIkR2dZJ8EuqP9qG9sddTvJx8wyJMruSUlkM2epCJa99cXzzwNlHjNys/o7eOip/Y/2h933EM4HzDIkytlszGF3YxsVK1FT848kzSQIjXfv2rOBMwcU6r5c3onmJ32EM4HrK4hV8pmYwo3ky1cUvh92seQbWoiapK2dv8ZXY3K9AZvJz6EvY4jebKF2ZUwZqdB3EK5v1cPnxNOsv7PyeWqP682YS2c54jqm+fQ8wkjHx7CTmRrP3kt7CfvHkbuWZSvNbrgxwnlktn8nq24/k1HzmLfP6+gP+EvuFzj2LLFTuUaQfoDHb3wU3/HBX4fhhT40dnTN3jPgPMfzG78e1brJ8+RPFku20oYt5RLysjKHk9/0o7nvluR8XGf+24FpnytxNB7ozZh7fch6YGh0JMC0tvewc2/R7dikCfLZVsJ45ZySRnZpOSe01cw5WslWd2D0fdGbcJaNpIXBX4RrQlmt/8e3YoTr2Q52WSb3kk4t5RLyqilQWqPXcjq2EbfG7UJ63LJ70P2daPc/nt0K47kKSt6cs3ZLghyU7kkkPyeaK3ybI/EslogZPS90UqrpP2eCsybLHXb79ErGOQpY3pzrNlWwpi9atRKqe9JR0+f5s9kk67I5L2RpVVSWw77fV9V1yR+P1Nu+j16CatrTODG2fhslZYOReX/+qupfVnUuKW6xkir3kRqWwVqMfu9MaMayq5rtYIb/57VqmsY5E3gxn8U2SotHYqKdYcg+sfjA/C+RsmdG+n5PU+rPS58T4z44ZTyrKpusqVWZmn2w9uJ3Pj37JgSyv7+fqxfvx5NTU0oLCzEK6+8gjFjxth5CWQi5liT1Te2wucDsh026am6sWpEXN/YKv0kwglSd7K1uubIkSPo7e3F22+/jerqamzatMnO05PJnNZaQLZc365zb/zLOWG5YdAHDCkw9qemVnUjaj8g2nPVKOW4Mvn68HY7W0fyH374ISorKwEAt912G86cOWPn6clkTmotoDYJbMc1ymrh/T6gZt4EzJtYljb67o72JW2Tl0it6saqenO1JmNqD2835Nnzma1BvqurC0VFRYP/HwgEEIvFEAwOXEZRUQjBoPs2FggE/CgtHZrry7CVcs9L7rwZS+682ZJz1J1uRu3hs2hpj2BUSRjVVRW4d8po4Wu3N1wWBr7//beL6I72IRJNCP6Hz2HY0JD0WDKy33Pd6WZpiiMex+D7k/pe1Z1uRvX//af0fNsbLgvfW7V682z+HaqlYzb84BvC96vudDM2Hj5nyvvrFF77e7Y1yBcVFeHatWuD/9/f3z8Y4AGgq8udOT83TtRky+p7Th2ZN7dHsPa/z+Da9R7hKLGlPSI8ztXr0bSvRaL9eO3PTbpa6CYS3bOeFIfofVJGv2pa2iPCn1WbC8nmd6K2i9XMMaXCY7/256bBAK/I9P11Cjf+PatNvNqak//Wt76F48ePAwBOnTqFiorcVRCQc4hy6WopCRGj+WKzJhEzSXEk5tTVyO7JqrkQ4XE1FkNxFavz2TqSr6qqQkNDAxYvXox4PI6NGzfaeXr6kpNyqLJcuixwqi3XF9V2FwZ8wgVJ2U4iam11B0BaV653g43rveK8vFVzIaLjam0awgor57M1yPv9frz00kt2npJSOK0ToGzELuuIaHS5PiBYqq9j1Ct6ECr5cT1b3ZUXh6Tvp95RbkdPn/R3Y+ZuU2rH1UpdcBWr87GtQZ5xWidAWcDrjw8Ei9RrTR3d6v1UYmTUK3sQDhs6kJvWGolrBTnZ6NcHpC2kcnqXRidVWJEYg3yecVoOVW2yb0XlWNQeu5BUZpg4ugWgu3eOkaAjexDWHj6LmT+epvpeaW3aAYhHv0B6gFc4Pb9t1acKMofmxOtPfvITHDlyBH192o2WyPmybftrNrVJxHkTyzCkIL2kVhndGp2c1UsWVJUKHtl7pSz71wp48yaW4Z5JI3VfD/PblA3NIL969Wr84x//wP3334/XXnsNly5dsuGyyCpOW6U6b2IZ1nzvFpQXh+DDQKBMnLBU++Rh1acSWXvgkiEDH3zNeA8bLl7V/Vrmtykbmuma8ePHY/Xq1fj888+xYcMGfP/738ftt9+Op556Ct/85jftuEYykRNzqGof99WqN2QrRoeHM89C1je2CjfIBoBrvX1J8wHZvId6H0Ql4SBTIZQVzb+Gd999F/v27cOFCxdw3333Yc2aNYjFYnjsscdQV1dnxzWSydyUQ1Wr3thy9LzwZ7JprLr1xCVEJfvdRfvig5Og2b6HatvtKcJBP6pnj8/4HESAjiBfV1eHJUuWYPr05BajP/vZzyy7KHKe1CoWrfpps6iNmmsONgl/Rs9GHTJaI2yt7+ut9hE9vII+YFgoiI5IzBGfsMgbNIN8bW2t8OtVVVWmXww5k6ikcO3+M1hTlf0mEnrIRs2yWvpsaI2w1SZBjaxBcGLajLyJJZSkSVjFEhXXb9u5mlYtwGe6b6qsvBHQXuJvdA2Cm9Jm5F4M8g7kpLYDgHZtvWyJv9WractVRt2ZLiAS7XPaHx84l1aKymlrEIgABnnHyUXbAa2HiiyFEQdw968b0B3tl05WWrlic0XlWGlePpvAKhthay3xZx8XciJbu1CSNqsW+Mjo2WVIVBeu6OjpkwZ4hVUj2XkTy1AiKZfMRWB12hoEIoAjecex+yO/njxyagrDKCXgmpWGSjzO8HAQQR8QS3jO5CqwcjKVnIhB3mHs/sivN9+uBC2jlIBrVhoq9TjtkRgK/D4ML/Sjs6cv54GVk6nkNAzyOaA2orW7davaQ0UUmI1IbNa1YMdJU7pfij55RPvjuKkwiKM/nWHo+ojyAYO8TmamGtRGtHZ/5Fd7qOjd3CJVMGHzaoVZaShWsBAZwyCvg5kVL3pz4EZTGJk+FNQeKi9IKlcAYHRJGC3tEfgEC5Ji8fQSRrPSUKxgITKG1TU6mFnxYvZIVE91jJZ5E8twYPl0vF89M6lVrlpL3Xd/8R28Xz0TsjYxqfczY9wNwtfJvi4jq/RRNhMhomQM8jqYGZjN7uduZcmlnpJAvfcja61rpOUu8FVr4tTSSWUzEQZ6omQM8jqYGZjNrqW2KketpICU/VaB9F7vgP77MfM6tTYTIaKv2JaT7+zsxKpVq9DV1YVoNIrnnnvONf3ozax4MXtiNdsctSifDyRvq6fstyq6Tr33Y3YunROwRPrYFuTffPNN3HHHHVi2bBkuXryI6upq7Nu3z67TZ8XswGxmLfWKyrF4+dDZpFWnBX6frgeQaEJZ1iYg20Zbeh+UeieROQFLpI9tQX7ZsmUoLCwEAPT19SEUctcfYy4XuWgFvtRNMvRummG0RDLbfjDKOWX3YaSKye71BERu5Ytns42OxO7du7Fz586kr23cuBGTJ09GW1sbHnvsMaxZswbTpk1Lek13dy+CQfH+mk4WCPjR16cdLOtON6P28Fm0tEcwqiSM6qoK3DtltObPrN1/BpGELenCBX5suO8buHfKaMza8jc0f7nBdKLRJWG8+4vvqB67Yt0hGPnlJx5T7z0bYfRejLyfmbz3qay4Z6fjPbtDgWCOSmFJkJdpamrCypUrsXr1asyaNSvt+21tnXZdiqm0uhMC6aNUYGDkmTqRmWrBjpPCtER5cQgHlk/HtNrj0kDtA1RTHrJji6Req557Nkp2Lz4A71fPzPi4mb73qay4Z6fjPbvDiBHF0u/ZVl1z/vx5PP3006itrRUGeK+TlTpuOXoeC3acxLTa41iw42RaCaDWBKNaDlqpm6852ITbBcdX6y6ZSFRVYwWzy0sVdnf2JHIS23LytbW16O3txYYNGwAARUVF2LZtm12nzzlZsO7o6Rvck1SUg5ZNMBaHAoZG4qLja3WXzGS0mw2r8uysxKF8ZluQz6eALqK1d6gitYpFtuFzd7QfHT3Gg5RaC4Vc70hlVd8eVuJQPmPvGoukBswZ427AO//6VFc1S+IIUxT4uqN9aI/E0n6u/MugpfUwkY1gndAm14prYCUO5TOueLWAqJ/MO//6FPdMGony4hB8GAjIenc1Su0t0yEI8MBA8NaTZ/f5kFfL/5VWCInvvZ1pKKJc4kg+Q4kj9VElYTw+Y0zSqFs00ddw8SoOLJ+edIxMRphq6Qc9uzj1x2H5vrFO44RPKUS5wCCfgdTg3NweSQqaeif6Ms1Ba6UfUvPs6+ub0toBW7nBtpZc5/6J8gmDfAZkI/X19QMtAYxM9IkmPl842KQa/Iw8HNT6wueiusTM3vxEpI05eYn6xlZp/bosOCppkBnjbjDcadJoX/h5E8uwonIsyopDaO3swdYTl6SvHa4z928H1qwT2YtBXkAr4KoFRyX3bnSiz2jw0/tQqG9sxbWe9IlavU3MzMaadSJ7MV0joLVFnygnnqi1sydtok/5ZCBLrxgNfnq2EVReFxP0ChhS4Fd96FiVN2fNOpG9OJIX0Aq4SkmesplGqtSAJRp11xxswqYjZ6U/o/V1PQ+F+sZWaYVN55erbEXM2FJQxuxNU4hIHUfyAsWhwGCrgUSJAVcZ1eopgZS19N1z+gqmfK1E+ulAFvzqG1vh80G4v6pyjUqglkm8l9RR+/XemK5PCZmwalUrEYkxyKeob2xFdzQ9IAd9SAu4SmDa3nAZLe0RacBSyzcrgVP5mdpjFwZXsxYG0j8qKME7tSQSSH4oaPWKv9LZgwU7TqatxFVbLWtW3pw160T2YZBPsfXEpaRdlhTDQkFpieKSO29WbU2q1rcmNXD2JARmZXNq5TzK9YmCt9+HpMldPQH5SmcP9py+ovk6BfPmRO7j2Zy8WgmkGmm3SEkrAT3U8s2JgVNPhY3s+uLx5DpzswMy8+ZE7uTJIJ/NxKEVPc3nTSzD7f81PO3rqYFTz2Sq3uvT2ytepiQcZK8XIg/wZLpGb3mhiBUdC+sbW/H/WrrSvn7PpJFpo2+t8kK91yea4LzeGxNOKKcKB/2onj2eQZ3IAzwZ5I3WnKdWl9wzaSQaLl41rfpDlkdvuHg16f/1BHCjLQ3UNspWjm/2/RKRc3gyyBtZcCPqpbL/n1cwLJT5W5P60NA76ao3gGdancLyRaL848kgbyTlIhplx+IYLGM02kBL9NCQ0WpYZgWWLxLlF09OvBrZJEJPqWEk1o/aYxd0nVurPj2RUqueTxt4EJG9PDmSB/SPWPXuvdoeiaG+sVXzmGoPjXLBuTL5pMB0CxHpZftI/sKFC5g6dSp6MtiE2gpGSg1lHSHrTjcP1uT7JP1syotDOLB8+uA+rIn0ttq1sqcMEXmTrUG+q6sLmzdvRmFhoZ2nVZWa2hkeCkhfKxql1ze2Yu3+M4OBV9ZuYMa4G7Bgx0ndk7AiW46eZy92IjLEtnRNPB7HunXrsHLlSqxYscKu0+qSmtq5+9cNmg3KFFtPXEJE0OvG/2UDsbLiUFp/GBGtxVb1ja3SGnf2YiciGUuC/O7du7Fz586kr40ePRrz58/HrbfeKv25oqIQgkH5SNqIutPNqD18Fi3tEYwqCaO6qgL3Thmt62dfWDAJa/efSQre4QI/Vs2ZgNLSoUmvVWszcPbluQCAWVv+phrgZcdOtL3hsvR7o0rCqj9rhUDAb/s5c433nB+8ds+WBPkHHngADzzwQNLXqqqqsGfPHuzZswdtbW149NFHsWvXrqTXdHWZMyIVbbS99r/P4Nr1Hl2TlNeu96DQ70Pky/8vCQdRPXs8Zo4pTWtEplaTr7y2pT2S9n1F+ZeTp6JjJ1I7xuMzxqj+rBVKS4fafs5c4z3nBzfe84gRxdLv2ZauOXz48OB/z549G2+88YZl58qmrYFoVWiPyih8ReVYbDx8LnnUn1KTL3sQKJOxesiOURIWd8ckIgI8WiefzT6iRvdanTexDBvu+4ZqTb4ZuyHJjlE9e7zuY8hk2rGTiJwvJ3Xyx44ds/T4slFvHMB3/8/fVZtvZfKAuHfKaMwcUyr9vhntBKxqSSBaoWukbp+InM2Ti6FWVI7FCwebIKhmRHskhpcPDeytKgpiVm00bUY7AStaEmST2iIi5/Nkuub0J+3CAK+I9sel6RezNpp2Swokm9QWETmfJ0fy+/6pvaWdLIiZkRZxUwrEqk8uROQMrg/yol4uolWnqYpDASzYcVIYyLNNi7gpBWLFJilE5ByuTtfIerlI2sck6ejpS/q5moNNuPvXDaakVdyUAjHSsZOI3MfVI3nZiDkc8CHSp2M4n6Kjp8+UtIrbUiDsMU/kXa4eyctGxpG+eNpoXqv52ODPmtDwy6zJWyKibLl6JC8bMft96d0g44CuTayBzNMqifMDxaEAQsEgOiIx9n0nopxxdZCXTRrq3ZlJJpO0SmpFTUdPH8JBP16cP4HBnYhyxtXpGtGk4T2TRsIvmXktCQc1NwjJNK1itB0CEZEdXD2SB5InDZXRtKyEsj0Sw/CUNMqMcTeg4eLVrFsFuKmihojyh+uDfCI9m2hblUaRzQ8MDwel9fhERFZzdbomld5RsxVpFFFFTYHfh2s9Me7JSkQ546kgb2TC1Ow0imh+YEiBH7GU1BHz9ERkJ0+la1ZUjkXNwSZdrxU9EEQtEoy2A058/bTa48LXMU9PRHbx1Eh+3sQylIS1n1uiChpZiwSt1Ep9YytmbfmbsNuk7JOFU1e+EpH3eCrIA0D17PHC3PjwUEC1N0smJZDKg6G5PSJ8MHDlKxHlmqfSNUDmrYIzKYHU6jZp1W5ORER6eSLIi3LpejfIVmTSVEzPg4HNv4gol1yfrsk0l54qk9QKc+5E5HS2Bfm+vj688sorWLx4Me6//3789a9/NeW4ZrUTyKSvOnPuROR0tqVr9u/fj1gshj/96U9obW1FfX29Kcc1s52A0dSK8trtDZfR0h5hzp2IHMe2IP/ee+/hlltuwfLlyxGPx7Fu3TpTjpvrDTrmTSzDkjtvxhdfXLflfERERvji8bjxLZQ07N69Gzt37kz62g033ICvf/3r2LhxIz744AP86le/wq5du5Je093di2BQe2OPRHWnm7F2/xlEognthgv82HDfN3DvlNGZ34QBgYAffX3ZtTd2G95zfuA9u0NBgTxuWhLkRX7+859j7ty5mDNnDgBgxowZaGhoSHpNW1tnRsfOdqVqtscrLR2adyN53nN+4D27w4gRxdLv2ZaumTp1Kt59913MmTMHH3/8MUaNGmXasc0sU0zd/EOp1lHOQ0TkJrZV1zz44IOIx+N48MEHsW7dOrz44ot2ndoQbv5BRF5i20i+sLAQr776ql2nyxg3/yAiL/HEilcz6anWSc3Zr5ozATPHlNp4lURE+rh+xavZtBY4iVbYrt1/hhuBEJEjMcin0Fr5KszZR5mzJyJnYrpGQK1ahzl7InITjuQNYlMyInITBnmDhDn7AjYlIyJnYrrGINFGIKyuISKnYpDPQGrO3o3LoIkoPzBdQ0TkYQzyREQexiBPRORhDPJERB7GIE9E5GEM8kREHsYgT0TkYZ6skzd7O0AiIrfyXJDn9n1ERF/xXLqG2/cREX3Fc0GerYCJiL7iuSDPVsBERF+xLch3dnbixz/+MR566CEsW7YMbW1tlpxHa/s+IqJ8YluQ37t3LyoqKvDWW29h/vz5+N3vfmfJebS27yMiyie2VddUVFTg4sWLAICuri4Eg9adWm37PiKifGJJpN29ezd27tyZ9LWamho0NDRg/vz5aG9vx65du6w4NRERJfDF4/G4HSf66U9/im9/+9tYvHgxPv74Y6xatQoHDhxIek13dy+CwYAdl2OqQMCPvr5+7Rd6CO85P/Ce3aGgQB43bUvXDB8+HMXFxQCAm266CdeuXUt7TVeXO8sc83FnKN5zfuA9u8OIEcXS79kW5J9++mk8//zzeOuttxCLxfDyyy/bdWoiorxlW5AvKyvD66+/btfpiIgINubkiYjIfp5b8UpERF9hkCci8jAGeSIiD2OQN9GFCxcwdepU9PS4sxTUiM7OTjz++ON45JFHsGjRInz00Ue5viTL9Pf3o6amBosWLcLSpUtx+fLlXF+S5aLRKFatWoWHHnoICxcuxNGjR3N9Sbb47LPPMGvWLFy4cCHXl2Iaz20akitdXV3YvHkzCgsLc30ptnjzzTdxxx13YNmyZbh48SKqq6uxb9++XF+WJY4cOYLe3l68/fbbOHXqFDZt2oRt27bl+rIsVVdXh9LSUrz22mv44osv8IMf/AB33313ri/LUtFoFDU1NQiHw7m+FFNxJG+CeDyOdevWYeXKlRgyZEiuL8cWy5Ytw+LFiwEAfX19CIW828r5ww8/RGVlJQDgtttuw5kzZ3J8RdabO3cunn76aQAD/74DAfetRDdq8+bNWLx4MUaOHJnrSzEVR/IGifryjB49GvPnz8ett96ao6uyluieN27ciMmTJ6OtrQ2rVq3CmjVrcnR11uvq6kJRUdHg/wcCAcRiMUub7OXasGHDAAzc+1NPPYVnnnkmtxdksb179+LGG29EZWUlduzYkevLMRXr5E1QVVWF8vJyAMCpU6cwefLkvGjA1tTUhJUrV2L16tWYNWtWri/HMq+++iqmTJmC+fPnAwBmzpyJ48eP5/iqrNfS0oInn3xyMC/vZQ8//DB8Ph98Ph8aGxsxduxYbNu2DSNGjMj1pWUvTqa666674pFIJNeXYblz587F58yZE29sbMz1pVju0KFD8WeffTYej8fjH330UfxHP/pRjq/Iem1tbfG5c+fG//73v+f6Umz3yCOPxM+fP5/ryzCNdz9vkqVqa2vR29uLDRs2AACKioo8OxlZVVWFhoYGLF68GPF4HBs3bsz1JVlu+/bt6OjowNatW7F161YAwOuvv+65Scl8wHQNEZGHsbqGiMjDGOSJiDyMQZ6IyMMY5ImIPIxBnojIwxjkiYg8jEGeiMjDGOSJVOzatQsrV64EADz77LN50a6CvIWLoYg0rFixAsOHD0dvby9++ctf5vpyiAxhkCfScOrUKSxatAh79+7FpEmTcn05RIYwyBOp6O3txSOPPIIf/vCH2LNnD/7whz/kzcYw5A3MyROp2LJlC77zne9g0aJFqKysRG1tba4vicgQjuSJiDyMI3kiIg9jkCci8jAGeSIiD2OQJyLyMAZ5IiIPY5AnIvIwBnkiIg9jkCci8rD/D0SB0X10lxJtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the data we have\n",
    "plt.scatter(X, y)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conditional probability of $y$ given $x$\n",
    "\n",
    "Using the linear model given above, compute the conditional distribution $y$ given $x$, $P(y \\mid x)$, as a function of $w$. Answer the following questions and __fill in the code below__.\n",
    "- __What family of distribution is the conditional probability?__\n",
    "- __What are its mean and variance?__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conditional(y, x, w):\n",
    "    \"\"\"\n",
    "    Arguments: y of shape (n, ), x of shape (n, ), and a scalar w.\n",
    "    Return: conditional probability of shape (n, )\n",
    "    \"\"\"\n",
    "    # TODO: Evaluate the conditional probability of y given x at some w.\n",
    "    ### start 1 ##\n",
    "    return np.exp(-0.5 * ((y - x.dot(w)) / sigma) ** 2) / (sigma * np.sqrt(2 * np.pi))\n",
    "    ### end 1 ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maximum likelihood (ML)\n",
    "\n",
    "We want to find the best weight $w$ that estimates the true unknown $w^*$. Given $n$ data samples $\\{(x_1, y_1),\\dots,(x_n, y_n)\\}$, we can compute the likelihood function:\n",
    "\n",
    "$$\\text{Likelihood}(w) := p(y_1,\\dots,y_n \\mid x_1,\\dots,x_n) = \\Pi_{i=1}^n ~p(y_i \\mid x_i)$$\n",
    "\n",
    "ML solution is the weight $w$ that, as the name suggests, maximizes the likelihood of the observed data. Run the following code to plot the likelihood and the log likelihood functions. Likelihood tends to get very small very quickly so often we compute the log conditionals and sum them instead.\n",
    "\n",
    "$$\\text{Log-Likelihood}(w) := \\log p(y_1,\\dots,y_n \\mid x_1,\\dots,x_n) = \\sum_{i=1}^n ~ \\log p(y_i \\mid x_i)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a7d93da5ad254c72a336cae728c96145",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=25, description='n', max=50, min=1), Output()), _dom_classes=('widget-in…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 200\n",
    "W = np.linspace(-5, 5, N)\n",
    "\n",
    "def compute_likelihood_1d(X, y, W, n):\n",
    "    likelihood = np.zeros(N)\n",
    "    for i in range(N):\n",
    "        likelihood[i] = conditional(y[:n], X[:n], W[i]).prod()\n",
    "    return likelihood\n",
    "\n",
    "\n",
    "def compute_loglikelihood_1d(X, y, W, n):\n",
    "    loglikelihood = np.zeros(N)\n",
    "    for i in range(N):\n",
    "        loglikelihood[i] = np.log(conditional(\n",
    "            y[:n], X[:n], W[i])).sum()\n",
    "    return loglikelihood\n",
    "\n",
    "\n",
    "def plot_likelihood_1d(n):\n",
    "    likelihood = compute_likelihood_1d(X, y, W, n)\n",
    "    loglikelihood = compute_loglikelihood_1d(X, y, W, n)\n",
    "    plt.figure(figsize=(10, 5))\n",
    "    \n",
    "    w_ml = W[likelihood.argmax()]\n",
    "    ml = likelihood.max()\n",
    "    \n",
    "    plt.subplot(121)\n",
    "    plt.plot(W, likelihood, color='blue')\n",
    "    plt.vlines(w_ml, 0, ml, color='black', linestyle='--', label=r'$w_{ML}$')\n",
    "    plt.xlabel('w')\n",
    "    plt.ylabel('Likelihood')\n",
    "    plt.xticks(np.arange(-5, 5, 1))\n",
    "    plt.legend()\n",
    "    \n",
    "    plt.subplot(122)\n",
    "    plt.plot(W, loglikelihood, color='blue')\n",
    "    plt.vlines(w_ml, loglikelihood.min(), loglikelihood.max(), color='black', \n",
    "               linestyle='--', label=r'$w_{ML}$')\n",
    "    plt.xlabel('w')\n",
    "    plt.ylabel('Log Likelihood')\n",
    "    plt.xticks(np.arange(-5, 5, 1))\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    \n",
    "    print('ML solution: w_ml = %.8f' % w_ml)\n",
    "\n",
    "interactive_plot = interactive(plot_likelihood_1d, n=(1, 50))\n",
    "output = interactive_plot.children[-1]\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__If your code is correct, with $n$ equal to 1, you should see a wide normal distribution with a peak density of around 0.4.__\n",
    "\n",
    "Answer the following questions:\n",
    "- __What is the ML solution? In other words, which value of $w$ maximizes the likelihood?__: $w_{ML} = 1$\n",
    "- __What do you notice when $n$ increases/decreases?__ (also see the y axis): Large $n$ makes distribution more peaky or concentrated at $w^*$, and the likelihood gets smaller as $n$ increases.\n",
    "- __Why does the shape of the log likelihood curve not seem to change with $n$?__: The log-likelihood curve just looks the same because of the log scale."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (c) Linear Regression (Posterior)\n",
    "\n",
    "Using the same model and data as in part (b), we want to take a different approach in finding the optimal $w^*$. Let's now assume that $w$ is a random variable drawn also from a Gaussian distribution with mean $\\mu_w$ and variance $\\sigma_w^2$, i.e. $w \\sim \\mathcal{N}(\\mu_w, \\sigma_w^2 I)$. This is an extra condition we impose on $w$, a _prior_ knowledge we have about $w$ before even seeing any data samples! So in this case, the _prior_ distribution of $w$ is $\\mathcal{N}(\\mu_w, \\sigma_w^2 I)$.\n",
    "\n",
    "The posterior distribution of $w$, given the same training data, is \n",
    "\n",
    "$$p(w \\mid y_1,\\dots,y_n, x_1,\\dots,x_n)$$\n",
    "\n",
    "It is the distribution of $w$ now that we see the training data. From Bayes' theorem, we can write the posterior in term of the likelihood and the prior.\n",
    "\n",
    "\\begin{align}\n",
    "p(w \\mid y_1,\\dots,y_n, x_1,\\dots,x_n) ~&=~ \\frac{p(y_1,\\dots,y_n \\mid x_1,\\dots,x_n,w) \\cdot p(w \\mid x_1,\\dots,x_n)}{p(y_1,\\dots,y_n \\mid x_1,\\dots,x_n)} \\\\\n",
    "&=~ \\frac{p(y_1,\\dots,y_n \\mid x_1,\\dots,x_n,w) \\cdot p(w)}{\\int_w p(y_1,\\dots,y_n \\mid x_1,\\dots,x_n,w) p(w) dw}\n",
    "\\end{align}\n",
    "\n",
    "We will derive the posterior analytically in the next section as well as in Problem 1 of the homework. For now, we will plot it and see what it looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "91e199421d404faa9c8156922c86d7e3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=13, description='n', max=25, min=1), IntSlider(value=0, description='mu_…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 200\n",
    "W = np.linspace(-5, 5, N)\n",
    "\n",
    "def plot_posterior_1d(n, mu_w, sigma_w):\n",
    "    likelihood = compute_likelihood_1d(X, y, W, n)\n",
    "    prior = np.exp(-0.5 * ((W - mu_w) / sigma_w) ** 2) / (sigma_w * np.sqrt(2 * np.pi))\n",
    "    py = simps(likelihood * prior, W)\n",
    "    posterior = likelihood * prior / py.prod()\n",
    "    w_map = W[posterior.argmax()]\n",
    "    \n",
    "    fig, ax1 = plt.subplots(figsize=(6, 6))\n",
    "    ax1.plot(W, posterior, color='blue', label='Posterior')\n",
    "    ax1.vlines(w_map, 0, posterior.max(),  color='black', \n",
    "               label=r'$w_{MAP}$', linestyle='--')\n",
    "    ax1.set_xlabel('w')\n",
    "    ax1.set_ylabel('Posterior')\n",
    "    ax1.set_xticks(np.arange(-5, 5, 1))\n",
    "    ax1.tick_params(axis='y', labelcolor='blue')\n",
    "    ax1.legend(loc=2)\n",
    "    \n",
    "    ax2 = ax1.twinx()\n",
    "    ax2.plot(W, prior, color='red', label='Prior')\n",
    "    ax2.set_ylabel('Prior')\n",
    "    ax2.tick_params(axis='y', labelcolor='red')\n",
    "    ax2.legend()\n",
    "    \n",
    "    fig.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "    print('MAP solution: w_map = %.8f' % w_map)\n",
    "\n",
    "interactive_plot = interactive(plot_posterior_1d, n=(1, 25), \n",
    "                               mu_w=(-5, 5), sigma_w=(0.1, 2, 0.2))\n",
    "output = interactive_plot.children[-1]\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interactive plot may have some delay. Give it some time :)\n",
    "\n",
    "Answer the following questions:\n",
    "- __Set mean of the prior $\\mu_w = 1$, what is the MAP solution?__: When $\\mu_w = 1$, there is no skew because $\\mu_w = w^* \\approx w_{ML}$.\n",
    "- __Set $\\mu_w = -1$ with a large $\\sigma_w$ and large $n$, what is the MAP solution?__: $w_{MAP}$ may move very slightly towards the mean of the prior $\\mu_w$, but the change is so small that it is difficult to see.\n",
    "- __With a small $\\sigma_w$ and small $n$, how does changing $\\mu_w$ affect MAP solution?__: With small $\\sigma_w$ and $n$, prior has a greater effect on the posterior. MAP solution is skewed towards the mean of the prior, even though the true $w^* = 1$.\n",
    "- __What is the effect of $n$ and $\\sigma_w$ on the posterior?__: Larger $n$ puts more emphasis on the likelihood term so it makes MAP solution closer to the ML solution or $w^*$. Intuitively, with more samples, we are more \"informed\" and more confident in $w_{ML}$, and relatively the prior knowledge becomes less important. When $\\sigma_w$ is large, we have a lot of uncertainty over our prior, which implies that the prior is \"uninformative,\" so once again the likelihood dominates the posterior.\n",
    "\n",
    "Well done! You have reached the end of the Jupyter notebook part. Go back to the worksheet and get a head start on the analytical parts!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
